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Introduction 

Optoacoustic  tomography  (OAT)  is  an  emerging,  hybrid  modality  that  offers  high  contrast  of  optical 
imaging  along  with  high  resolution  of  ultrasound  imaging  [1].  While  several  groups  have  looked  at  the 
application  of  OAT  in  detection  of  breast  cancer  [1]  and  optoacoustic  molecular  imaging  (OMI)  [3]  no 
one  has  yet  looked  using  OMI  for  breast  cancer  using  proteases  as  target.  Elevated  levels  of  protease  are 
found  in  several  cancers  including  breast  cancer  [2] .  This  proposal  will  aid  in  the  design  and 
development  of  an  optoacoustic  molecular  imaging  system  specifically  targeted  to  detect  elevated  levels 
of  protease.  We  are  specifically  implementing  software  that  will  simulate  the  following:  ideal 
optoacoustic  signal,  effect  of  finite  transducer  width  and  temporal  response,  effect  of  non-uniform 
illumination,  effect  of  variable  speed  of  sound,  and  effect  of  ultrasonic  attenuation.  In  addition,  we  will 
compare  the  various  image  reconstruction  algorithms,  validate  the  simulations  using  experimental  data 
and  refine  the  optoacoustic  device  design. 

This  technique  has  the  potential  to  detect  breast  cancer  earlier  than  mammography,  especially  in 
young  women  with  radiographically  dense  breasts.  It  also  has  the  potential  to  allow  for  molecular 
characterization  of  tumors  in  women  of  all  ages.  This  technique  has  also  the  added  advantage  of  being 
non-invasive  and  involves  no  ionizing  radiation.  It  is  also  relatively  inexpensive  with  costs  comparable  to 
that  of  ultrasound 

This  report  summarized  the  progress  made  on  this  pre-doctoral  traineeship  project  for  breast 
cancer  by  the  recipient  during  the  past  year. 
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Body 


Training: 

The  program  requires  22  courses  over  all.  However,  the  applicant  needs  to  take  only  15  courses  since  she  has 
a  Masters  degree  in  Physics  and  had  previous  teaching  experience.  The  applicant  has  taken  14  of  the  required 
15  courses.  The  courses  taken  in  the  past  year  include  cancer  biology. 

Research  Accompiishments 

The  main  goals  of  this  project  are: 

(1)  developing  and  implementing  simulations  of  optoacoustic  physics  and  signal  detection  with  a  handheld 
optoacoustic  probe,  with  the  aim  of  performing  optoacoustic  molecular  imaging. 

(2)  validating  the  simulation  software  by  comparing  it  with  actual  optoacoustic  signal  generation  and 
detection  in  the  laboratory. 

(3)  comparing  image-reconstruction  software  for  optoacoustic  image  with  a  handheld 
probe. 

(4)  using  these  tools  to  refine  the  design  of  the  optoacoustic  image  device. 

We  performed  the  following  research  during  the  past  year: 

1.  Development  of  ideal  optoacoustic  signal  software 

We  have  developed  the  ideal  optoacoustic  signal  software  for  the  specific  2D  geometry  that  we  are 
considering.  In  this  geometry  the  transducer  is  scanned  along  a  line  and  the  inherently  3D  problem  is 
reduced  to  2D  because  of  transducer’s  directivity.  We  have  also  included  the  finite  length  of  the 
stimulating  pulse  and  the  finite  temporal  response  of  the  transducer  in  to  our  simulation  software.  Both 
these  effects  are  assumed  to  be  ideal  delta  functions  in  most  reconstruction  algorithms. 

2.  Studying  the  effect  of  variable  speed  of  sound  and  ultrasonic  attenuation 

There  are  several  physical  factors  that  impact  image  quality  in  OAT.  These  include  limited  data,  noise, 
non-uniform  illumination,  variable  speed  of  sound  and  ultrasonic  attenuation.  There  are  several 
investigators  looking  at  the  problems  of  limited  data  and  non-uniform  illumination.  We  investigated  two 
of  these  factors:  variable  speed  of  sound  and  ultrasonic  attenuation. 

The  majority  of  the  image  reconstruction  algorithms  in  optoacoustic  tomography  (OAT)  so  far  have 
assumed  a  uniform,  non-dispersive  acoustic  medium.  The  effect  of  frequency-dependent  attenuation  on 
acoustic  waves  can  be  significant  since  OAT  uses  broadband  detection.  Reconstructed  images  may 
exhibit  distortion  and  artifacts  if  these  effects  are  not  taken  into  account.  Previous  work  on  dispersive 
acoustic  media  done  by  our  group  [4]  focused  on  incorporating  the  frequency  dependent  attenuation 
effects  into  the  OAT  model.  We  explored  this  effect  further  and  studied  the  conditioning  of  the  resulting 
inverse  problem.  We  have  derived  an  inversion  formula  for  the  optical  absorption  function  in  the 
presence  of  attenuation  using  singular- value  decomposition  (SVD).  This  expression  relates  the  measured, 
attenuated  pressure  via  an  integral  operator  to  the  optical  absorption  coefficient  of  the  medium  [10].  We 
also  derived  another  expression  for  attenuation  correction  when  a  source  is  placed  in  a  medium  with 
different  attenuation  from  that  of  the  source.  This  is  a  common  scenario  in  practice  where  a  tumor  may 
have  a  different  attenuation  coefficient  from  that  of  the  surrounding  healthy  tissue. 

Speed  of  sound  can  vary  in  the  breast  tissue  by  as  much  as  10%.  The  variable  speed  of  sound  affects  the 
travel  times  of  the  sound  waves  that  reach  the  transducer.  Incorrectly  computed  travel-times  can  have  an 
effect  on  image  quality  and  accuracy.  Previous  research  in  this  area  has  only  looked  at  first-order  travel¬ 
time  perturbations[5,6]  and  has  not  incorporated  the  effects  of  ray  bending.  We  systematically  explored 
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this  effect  and  developed  algorithms  for  incorporating  it  into  our  imaging  model.  We  have  derived 
improvements  to  the  zeroth-order  geometrical  acoustics  (GA)  model  that  incorporates  second-order 
perturbations  into  travel  times  of  sound  waves  [11].  For  homogeneous  speed  of  sound,  the  travel-time 
surfaces  of  equal  times  (isochrons)  are  spheres  in  3D  and  circles  in  2D.  With  variable  speed  of  sound 
theses  surfaces  are  distorted  based  on  the  variation  in  speed.  Perceivable  differences  have  been  detected 
in  the  isochronous  surfaces  in  the  inhomogeneous  case  with  speed  of  sound  variation  of  10%  or  more. 


3.  Comparison  of  optoacoustic  image  reconstruction  algorithms 

There  are  several  algorithms  that  exist  in  OAT  for  the  2D  geometry  that  we  are  considering.  However, 
these  algorithms  have  not  been  evaluated  for  noise  and  resolution  properties.  It  is  important  to  know 
which  algorithm  may  be  best  suited  for  optoacoustic  molecular  imaging  of  our  molecular  probe.  We 
compared  three  image-reconstruction  algorithms  for  optoacoustic  imaging  with  a  handheld  probe.  Two  of 
these  are  based  on  previously  published  analytical  algorithms  that  assume  fairly  idealized  data  and  had 
already  been  implemented  by  us  in  the  past  year.  The  first  one  is  the  Fourier-based  algorithm  [7]  that  is 
exact  in  3D  but  not  in  the  2D  geometry  that  we  are  considering.  The  second  one  is  the  synthetic  aperture 
algorithm  [8]  that  is  approximate  but  is  the  most  commonly  used  by  researchers  in  the  field.  The  third 
algorithm  is  based  on  Norton's  approach  for  reflection-mode  tomography  [9]  and  has  been  applied  by  us 
to  OAT  for  the  first  time. 

We  found  that  the  2D  Fourier-based  algorithm  offers  better  resolution  and  local  modulation  transfer 
function  (LMTF)  in  the  depth  direction  while  Norton's  algorithm  offers  the  best  lateral  resolution  [12]. 
However,  we  found  that  in  reconstructions  of  phantoms  the  images  produced  by  Norton-based  algorithm 
looked  the  sharpest  and  more  uniform.  The  LNEQ  data  suggests  that  Norton-based  algorithm  has  the  best 
signal  detectability. 

We  will  also  implement  a  new  algorithm  that  we  have  derived  using  penalized-likelihood  iterative 
algorithms,  which  have  not  previously  been  applied  in  the  context  of  OAT  but  which  have  the  potential 
to  model  more  of  the  physical  degradations  inherent  in  the  modality  and  thus  to  lead  to  higher  quality 
images.  These  results  will  be  validated  in  the  lab  in  the  coming  year  when  the  molecular  probe  becomes 
available  ad  the  imaging  system  is  ready. 
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Key  Research  Accomplishments 

1.  We  have  developed  software  that  simulates  ideal  optoacoustic  signal.  We  have  also  incorporated  the 
finite  transducer  and  finite  laser  pulse  length  effects  into  this  simulation  software. 

2.  We  have  studied  the  effect  of  variable  speed  of  sound  and  have  derived  a  higher-order  geometrical 
acoustics  expression  that  addresses  this  effect  in  image  reconstruction. 

3.  We  have  derived  a  new  way  of  incorporating  the  effect  of  ultrasonic  attenuation  in  optoacoustic 
image  reconstruction. 

4.  We  have  performed  a  systematic  comparison  of  several  algorithms  for  2D  image  reconstruction  in  terms  of 
resolution,  contrast,  noise  and  signal  detectability  properties.  We  found  that  an  algorithm  that  was  derived  b 
us  based  on  Norton’s  approach  had  the  best  contrast,  resolution  and  signal  detectability  properties. 
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Conclusions 


The  recipient  of  the  pre-doctoral  training  award  has  one  more  required  course  left  towards 
completing  the  PhD.  requirements.  She  also  needs  to  take  a  course  on  molecular  biology  to  expand  her 
knowledge  of  the  life  sciences. 

During  the  first  year  of  the  award,  we  focused  on  implementing  the  foundation  of  the  simulation 
software  that  would  aid  in  the  design  and  development  of  the  optoacoustic  molecular  imaging  system  being 
developed  by  our  group.  We  also  systematically  studied  two  very  important  factors  that  impact  image  quality 
in  OAT:  variable  speed  of  sound  and  ultrasonic  attenuation.  We  derived  methods  to  incorporate  these  into 
image  reconstruction  and  developed  algorithms  that  can  correct  for  these  physical  factors.  We  also 
systematically  evaluated  three  image  reconstruction  algorithms  in  OAT  in  a  homogeneous,  non-dispersive 
medium  in  2D.  We  found  that  the  new  algorithm  that  was  applied  by  us  to  OAT  based  on  Norton’s  approach 
offered  the  best  contrast,  resolution  and  signal  detectability. 

Our  future  research  will  focus  on  validation  of  our  simulations,  recommendations  for  the  design  of 
the  OMI  system  and  its  refinement. 
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Reportable  Outcomes 


Peer  Reviewed  Journals 

"Implementation  and  comparison  of  reconstruction  algorithms  for  2D  optoacoustic  tomography  using  a 
linear  array",  D.Modgil  and  P.J.La  Riviere,  accepted  for  publication  in  the  Journal  of  Biomedical  Optics, 
2009 

Conference  Proceeding  Articles 

1.  Dimple  Modgil,  Mark  A.  Anastasio,  Kun  Wang  and  Patrick  J.  La 
Riviere,  "Image  reconstruction  in  photoacoustic  tomography  with  variable 
speed  of  sound  using  a  higher  order  geometrical  acoustics  approximation" 

(Proceedings  Paper)  Proceedings  Vol.  7177,  Photons  Plus  Ultrasound:  Imaging 
and  Sensing  2009,  Alexander  A.  Oraevsky;  Lihong  V.  Wang,  Editors,  7 1771  A, 

February  2009 

2.  Dimple  Modgil,  Mark  A.  Anastasio,  Patrick  J.  La  Riviere,  "Photoacoustic 
image  reconstruction  in  an  attenuating  medium  using  singular  value 
decomposition"  (Proceedings  Paper) 

Proceedings  Vol.  7177,  Photons  Plus  Ultrasound:  Imaging  and  Sensing  2009, 

Alexander  A.  Oraevsky;  Lihong  V.  Wang,  Editors,  7 177  IB,  February  2009 

3.  "Photoacoustic  image  reconstruction  in  an  attenuating  medium  using 
singular-value  decomposition",  Modgil,  Dimple  La  Riviere,  Patrick  J. 

Nuclear  Science  Symposium  C  Conference  Record,  2008.  NSS  '08.  IEEE  Pub.  Oct. 

2008,  pp  4489-4493 


Honors  and  Awards 

Biological  Sciences  Division  Travel  award.  University  of  Chicago,  Spring  2009  for  SPIE  Photonics  West 
meeting,  2009 
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Implementation  and  comparison  of  reconstruction  algorithms 
for  2D  optoacoustic  tomography  using  a  linear  array 

Dimple  Modgil  and  Patrick  J.  La  Riviere 
Department  of  Radiology,  The  University  of  Chicago 


ABSTRACT 

The  goal  of  this  paper  is  to  compare  and  contrast  various  image  reconstruction  algorithms  for  optoacoustic  tomography 
(OAT)  assuming  a  finite  linear  aperture  of  the  kind  that  arises  when  using  a  linear-array  transducer.  Because  such  trans¬ 
ducers  generally  have  tall,  narrow  elements,  they  are  essentially  insensitive  to  out  of  plane  acoustic  waves,  and  the  usually 
3D  OAT  problem  reduces  to  a  2D  problem.  Algorithms  developed  for  the  3D  problem  may  not  perform  optimally  in 
2D.  We  have  implemented  and  evaluated  a  number  of  previously  described  OAT  algorithms,  including  an  exact  (in  3D) 
Fourier-based  algorithm  and  a  synthetic  aperture  based  algorithm.  We  have  also  implemented  a  2D  algorithm  developed 
by  Norton  for  reflection  mode  tomography  that  has  not,  to  the  best  of  our  knowledge,  been  applied  to  OAT  before.  Our 
simulation  studies  of  resolution,  contrast,  noise  properties  and  signal  detectability  measures  suggest  that  Norton’s  approach 
based  algorithm  has  the  best  contrast,  resolution  and  signal  detectability. 

Keywords:  Optoacoustic  tomography,  photoacoustic  tomography,  thermoacoustic  tomography,  image  reconstruction 


1,  INTRODUCTION 

Optoacoustic  imaging  is  a  hybrid  imaging  technique  that  has  attracted  a  lot  of  attention  in  recent  years  [1,2].  It  is  based  on 
the  photoacoustic/optoacoustic  effect,  which  refers  to  acoustic  wave  generation  upon  absorption  of  pulsed  optical  energy 
by  a  medium.  A  slight  rise  in  temperature  of  the  medium  due  to  the  absorption  of  the  incident  electromagnetic  wave  results 
in  thermoelastic  expansion.  This  thermoelastic  expansion  and  subsequent  contraction  leads  to  the  generation  of  acoustic 
waves.  Under  the  constraints  of  thermal  and  stress  confinement,  this  thermal  expansion  leads  to  a  rise  in  pressure,  p{r,t), 
that  satisfies  the  three-dimensional  inhomogeneous  wave  equation  [3]: 


d^p{r,t) 


(3  d 


(1) 


where  the  heating  function,  is  the  thermal  energy  deposited  by  the  electromagnetic  radiation  per  unit  time  per  unit 

volume,  (3  is  the  isobaric  volume  expansion  coefficient,  and  Cp  is  the  specific  heat  of  the  medium.  The  heating  function 
can  be  expressed  as  the  product  of  a  spatially  varying  absorbed  optical  energy  in  the  medium,  A(r),  and  a  time  dependent 
optical  illumination  function,  I{t).  The  absorbed  optical  energy  in  the  medium  is  the  product  of  the  optical  absorption 
function  and  the  optical  fluence. 


Optoacoustic  tomography  is  inherently  a  three-dimensional  inverse  problem.  The  sound  waves  generated  by  a  3D 
distribution  of  optoacoustic  sources  are  spherical  waves  radiated  into  the  volume  surrounding  the  sources.  These  3D 
optoacoustic  signals  can  be  detected  using  isotropic  ultrasound  detectors  arrayed  on  a  2D  measurement  aperture,  and  a  3D 
image  can  be  reconstructed  using  these  signals.  However,  detection  of  these  signals  using  a  ID  linear  array  of  transducers 
and  reconstruction  of  a  2D  image  slice  is  sometimes  more  practical  and  cost-effective  especially  in  a  clinical  setting.  The 
problem  can  be  reduced  to  2D  by  making  one  of  the  following  assumptions,  using  the  terminology  in  the  paper  by  Kostli 
et  al.  [4]: 


1 .  two-dimensional  source  distribution 

2.  two-dimensional  source  directivity 

3.  two-dimensional  detector  directivity 
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2D  source  distribution  implies  that  the  source  is  truly  two-dimensional  and  it  lies  in  the  detection  plane.  This  is  not  a 
very  realistic  scenario  for  biomedical  applications  since  there  are  generally  3D  sources  present  in  the  human  body.  The 
second  assumption  implies  that  even  though  the  source  is  3D,  it  is  constant  in  the  third  direction.  This  kind  of  source  will 
be  highly  directional  and  the  signals  received  by  the  detector  will  only  be  from  sources  in  the  detection  plane.  The  third 
assumption  is  relevant  for  highly  directional  detectors  that  are  insensitive  to  signals  coming  from  outside  the  detection 
plane.  Thus  a  2D  cross-sectional  slice  of  the  unknown  source  is  reconstructed  and  the  image  reconstruction  problem  is 
reduced  to  2D.  All  three  assumptions  imply  that  detected  acoustic  signals  are  only  from  in-plane  sources.  However,  these 
assumptions  are  not  exactly  equivalent.  The  first  two  assumptions  impose  constraints  on  the  source  geometry  that  may  not 
exist.  This  limitation  can  affect  the  accuracy  of  the  resulting  2D  reconstructed  image  which  may  be  especially  important 
for  quantitative  optoacoustic  imaging.  In  this  paper,  we  consider  the  scenario  based  on  the  third  assumption  where  the 
problem  is  reduced  to  2D  due  to  detector’s  directivity.  This  is  achieved  by  using  a  linear  array  of  anisotropic  transducers 
which  have  tall,  narrow  elements  that  are  essentially  insensitive  to  out  of  plane  acoustic  waves. 

There  are  several  algorithms  that  have  been  proposed  for  3D  image  reconstruction  in  OAT  for  measurements  over 
a  planar  aperture.  These  include  the  Fourier-based  algorithm  [4—6],  synthetic  aperture  (SA)  algorithm  [7,  8]  synthetic 
aperture  plus  coherent  weighting  algorithm  [9],  and  universal  backprojection  algorithm  [10].  Some  of  these  algorithms 
have  also  been  applied  to  image  reconstruction  in  2D  OAT.  The  Fourier  based  algorithm  is  theoretically  exact  in  3D  for 
continuously  sampled  data  on  an  infinite  measurement  aperture  but  not  necessarily  in  2D.  This  algorithm  implicitly  uses 
the  second  assumption  above  when  applied  in  2D,  namely  that  the  object  is  does  not  vary  in  the  third  direction.  Synthetic 
aperture  and  coherent  weighting  algorithms  [7-9],  on  the  other  hand,  are  approximate  reconstruction  algorithms.  There 
exists  an  algorithm  in  reflection  mode  tomography  that  was  proposed  by  Norton  [11],  which  is  theoretically  exact  in  2D 
for  continuously  sampled  data  on  an  infinite  measurement  aperture.  We  have  applied  this  algorithm  to  OAT  for  the  first 
time  (to  the  best  of  our  knowledge).  We  have  implemented  this  algorithm  for  a  planar  geometry  using  a  linear  transducer 
array  with  tall,  narrow  elements. 

Of  course,  no  algorithm  is  theoretically  exact  for  sampled  data  acquired  on  a  finite  interval,  so  this  paper  compares 
these  algorithms  for  that  practically  relevant  regime  by  examining  image  contrast,  resolution  and  noise  properties.  The 
studies  are  performed  using  simulated  optoacoustic  pressure  data.  We  go  beyond  the  standard  image  quality  metrics  by 
computing  noise  texture  measures  like  local  noise  power  spectra  (LNPS)  and  resolution  measures  like  local  modulation 
transfer  function  (LMTF).  These  noise  and  resolution  measures  are  used  to  obtain  the  local  noise  equivalent  quanta  (LNEQ) 
metric  that  is  known  to  predict  signal  detectability  under  certain  conditions  [12]. 


2.  METHODS 


The  optoacoustic  pressure  signals,  p(r,  t),  for  an  impulse  optical  illumination,  are  related  to  the  absorbed  optical  energy 
A(r)  [13]  as 


p{r,t)  =  T] 
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where  rj  =  Eq.  (2)  states  that  the  time  integral  of  acoustic  pressure  at  a  point  r  and  time  t  is  given  by  the  integral  of  the 
absorbed  optical  energy  function  over  a  spherical  surface  of  radius  |r  —  r'|  =  cf  centered  at  r.  A  simple  but  inexact  way  to 
reconstruct  A(r)  is  to  spatially  resolve  the  optoacoustic  waves  by  using  the  speed  of  sound  and  to  backproject  them  over 
hemispheres.  In  2D  geometry,  this  reduces  to  spatially  resolving  the  optoacoustic  waves  by  using  the  speed  of  sound  and 
backprojecting  over  semicircles  to  obtain  a  two-dimensional  slice  of  A(r).  This  is  the  method  followed  by  the  synthetic 
aperture  algorithm. 


Let  us  consider  a  line  of  transducers  along  the  x  axis  (i.e.  at  y  =  0,  z  =  0).  Let  A{x,  z)  represent  an  effective  2D  slice 
of  the  absorbed  optical  energy  function  in  the  half-plane  (y  =  0,  z  >  0)  .  The  pressure  signals  in  2D  reduce  to  (as  derived 
in  appendix  A) 


p{x,z  =  0,t) 
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Define  g{x,  t)  as: 


g{x,t)  =  —  [  p{x,t')dt' 

VC  J 

=  11  (ci  -  .  (4) 

We  will  consider  3  algorithms  for  our  study:  an  approach  based  on  Norton’s  algorithm,  the  Fourier-based  approach, 
and  the  synthetic  aperture  algorithm.  We  will  not  consider  the  synthetic  aperture  plus  coherent  weighting  algorithm  since 
it  is  non-linear  due  to  the  presence  of  the  coherence  factor,  and  nonlinear  algorithms  cannot  be  meaningfully  characterized 
using  the  LMTF,  LNPS  and  LNEQ  functions. 

2.1.  Application  of  Norton’s  algorithm  for  reflection  mode  tomography  to  OAT 

In  this  section  we  derive  an  exact  expression  for  optoacoustic  image  reconstruction  in  2D  closely  following  derivation  of 
an  algorithm  by  Norton  for  reflection  mode  tomography  [11]. 

Letting  r  =  ct  and  using  the  identity,  5{r  —  a)  =  2rS{r'^  —  a^),  one  can  write  equation(4)  as: 


g{x,t)  =  2r  J  J  dx'dz'A{x',z')S{r^  —  {x  —  x')‘^  —  z'‘^). 

(5) 

Defining  new  variables,  p  = 

r^,  C  =  z'^,  and  substituting  in  equation  (5)  yields 

g{x,p)  =  y/pj  JdxdC  ^  5{p  (x  x )  (). 

(6) 

Finally,  setting 
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and 
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,  equation  (6  )  becomes 

g'ix,  P)  =  J  dx'dCA'{x,  C)6{p  -C  -{x-  x'f), 

(9) 

which  can  be  written  as  a  2D  convolution. 


g'{x,p)  =  A'{x,p)**d{p-x‘^).  (10) 

This  convolution  relation  can  in  principle  be  solved  for  A'  by  taking  2D  Fourier  transform  on  both  sides  of  equation 
(10)  with  respect  to  x  and  p.  The  2D  convolution  in  Fourier  space  becomes  multiplication  of  the  2D  Fourier  transforms 
due  to  the  convolution-multiplication  theorem.  This  can  then  be  explicitly  solved  for  the  Fourier  transform  of  A'  and  on 
taking  the  inverse  2D  Fourier  transform  one  obtains  A'{x,  p).  The  transformation  of  g{x,  t)  into  g'{x,  p)  by  equation  (7) 
can  be  regarded  as  a  ID  geometric  distortion  of  the  cc  —  r  plane  in  the  r  direction.  The  double  convolution  relation,  equation 
(10)  implies  that  the  time-integrated  pressure  signal  is  equivalent  to  a  linear,  space-variant  mapping  of  the  absorbed  optical 
energy  from  points  in  the  x  —  z  object  plane  to  hyperbolas  in  the  x  —  r  measurement  plane  (please  see  reference  [11], 
section  II  for  more  details). 

Another  approach  that  is  more  direct  to  solving  equation  (10)  is  to  seek  a  solution  such  that: 

A'{x,p)  =  g'{x,p)**R{x,p)  (11) 

where  we  need  to  determine  R{x,  p).  This  can  be  done  using  the  method  outlined  in  Norton’s  paper  [11]  and  this  is  also 
derived  in  appendix  B  and  the  final  result  is: 


g{x' ,  r)Ri  [vc{z‘^  —  +  (x  —  x')'^}]drdx' . 


(12) 


A{x,  z)  =  2zv‘l 


where  Ri{u)  =  4sinc(2t6)  —  2sinc^(M)  and  iZc  is  the  cut-off  frequency  that  dictates  the  band-limit  of  the  measured  signals. 

The  above  relation  is  an  exact  equation  relating  the  absorbed  optical  energy  in  the  medium  to  a  filtered  back-projection 
of  time-integrated  pressure  signals. 

_  1, 

In  the  case  when  z  >>  Vc  ^ ,  this  can  be  approximated  as: 


A{x,z)  =  zi'i  J  J  -I-  (a:  —  —  r}]drdx' . 

Defining  G{x' ,  r)  =  *  Ri['s/i'^r\  and  substituting  in  Eq.  (13)  yields: 


(13) 


A{x,  z)  =  zvc  J  G  ^x',  \/z^  +  (x  —  x')^'j  dx' .  (14) 

Equation  (13)  is  similar  to  the  filtered  backprojection  (EBP)  expression  used  for  image  reconstruction  in  computer 

tomography  (CT).  The  function  i?i  is  the  same  as  the  Eourier  transform  of  the  truncated  ramp  filter  used  in  the  EBP 

expression  [14].  Eq.  (14)  can  be  seen  as  backprojection  of  the  filtered  function  G.  So  this  algorithm  is  equivalent  to 

_  1 

ID  filtration  at  each  transducer  position  x'  followed  by  backprojection  on  circular  arcs.  Here,  Vc  ^  can  be  regarded  as  a 
measure  of  the  resolution  of  signal  g{x,  r)  in  the  temporal  direction  since  Vc  is  the  bandwidth  of  the  function  g' {x,  p)  with 
respect  to  the  square  of  the  temporal  variable  r.  We  found  that  the  exact  Norton’s  algorithm  was  extremely  sensitive  to  the 
choice  of  cut-off  frequency  and  did  not  give  us  good  results.  Hence,  we  implemented  the  approximate  equation  (13)  as 
Norton-based  algorithm. 

So  A{x,  z)  can  be  obtained  via  the  following  steps  in  the  approximate  Norton-based  algorithm: 

1.  Convolve  the  2D  time-integrated  pressure  signal  g  with  ID  filter  Ri  with  respect  to  the  temporal  variable  r  . 

2.  Map  the  result  onto  a  circular  grid  for  a  given  {x,z). 

3.  Sum  the  resulting  expression  over  all  the  transducers. 

4.  Multiply  the  result  by  the  distance  from  the  transducer  axis,  z  and  other  constant  factors. 

2.2.  Fourier-based  algorithm 

This  algorithm  has  been  derived  by  Kostli  et  al.  [4]  and  L.V.Wang  et  al  [5].  It  relates  the  Eourier  transform  of  the  absorbed 
optical  energy  function  to  the  Eourier  transform  of  the  measured  optoacoustic  pressures.  The  relation  in  2D  is  given  by  [4]: 

.4 

where 

P{kx,uj)  =  J  J p{x,t)  Gxp{—iko!:x)cos{u}t)dxdt, 

and  w  =  Cy/k'^  +  k^.  Note  that  our  notation  is  different  from  that  in  reference  [4]. 

Here,  A{x,  z)  can  be  obtained  via  the  following  steps: 

1.  Take  the  real  part  Eourier  transform  of  pressure,  p{x,  t)  with  respect  to  time. 

2.  Take  the  Eourier  transform  of  the  result  with  respect  to  x.  This  gives  us  P{kx,oj)- 

3.  Scale  P{kx,uj)  via  equation  (15)  to  obtain  A  [kx,  kz  =  \/(^)^  — 


(15) 

(16) 


4.  Map  A  (k,kz  =  \/(f)^  —  k^)  to  a  Cartesian  grid  via  interpolation  to  obtain  A{kx,kz).  We  employed  bilinear 
interpolation.  Note  that  only  values  for  which  k^  <  (-)^are  used  in  the  mapping.  Higher  values  would  produce 
imaginary  values  of  kz  (physically  they  correspond  to  rapidly  decaying  evanescent  wave  components). 

5.  Inverse  Fourier  transform  A{kx,  kz)  to  obtain  A{x,  z). 


2.3.  Synthetic  aperture  based  algorithm 

This  algorithm  relates  the  signal  intensity  at  each  image  point  A(x,  z)  to  the  sum  of  signals  from  the  transducer  at  different 
positions  delayed  with  the  the  time  it  takes  the  signal  to  travel  from  the  transducer  position  to  the  point  [8]. 
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So  the  steps  involved  for  obtaining  A{x,  z)  are: 

1 .  For  each  point  (x,  z)  in  the  source,  sum  over  time-integrated  pressure  samples  corresponding  to  transducer  positions 
x'  at  different  times  t  such  that  y/{x  —  x')^  +  z^  =  cl. 

2.4.  Details  of  the  simulation 

All  the  simulations  were  performed  using  128  transducer  positions  spaced  0.1  mm  apart  and  128  time  samples  of  width 
67  ns.  Time-integrated  pressure  data  was  simulated  by  integrating  the  optoacoustic  source  over  circles  of  radii  ct  with  the 
transducer  placed  at  the  center  of  these  circles.  Simulated  2D  pressure  data  was  used  for  the  Fourier  based  algorithm  while 
time-integrated  simulated  pressure  data  was  used  for  the  synthetic  aperture  and  Norton  based  algorithms.  In  order  to  focus 
simply  on  the  acquisition  geometry  and  the  inherent  differences  between  the  algorithms,  we  did  not  simulate  a  low  pass  or 
bandpass  transducer  response. 

Simulated  pressure  data  was  generated  for  two  different  phantoms  -  a  circular  phantom  and  a  phantom  consisting  of  a 
line  of  rectangles.  An  exact  analytical  expression  was  used  to  generate  the  time-integrated  pressure  data  for  the  circular 
and  point  source  phantoms  A  discrete  numerical  method  was  used  to  generate  the  time-integrated  pressure  data  for  the 
line  phantom.  This  method  is  based  on  an  implementation  of  a  variation  of  Siddon’s  algorithm  [15]  for  computing  the 
intersection  lengths  of  an  arc  specified  by  the  coordinates  of  a  source  and  receiver  with  a  pixel.  The  circular  phantom  was 
of  radius  1  mm  with  its  center  placed  at  a  distance  of  2  mm  from  the  transducer  axis.  The  line  of  rectangles  was  placed  at  a 
distance  of  2  mm  from  the  transducer  axis  with  each  rectangle  being  0.5  mm  x  0.3  mm  wide.  The  images  were  constructed 
on  a  128  x  128  grid. 

To  study  the  resolution,  simulated  pressure  data  was  generated  for  a  point  source  of  size  0.1  mm  (same  as  pixel  width) 
placed  at  a  distance  of  1 .0  mm  from  the  transducer  axis.  A  zoomed-in  image  of  a  point  source  of  was  reconstructed  using  the 
3  algorithms  with  a  zoom  factor  of  10.  The  images  were  reconstructed  on  a  64x64  grid.  We  used  these  point  source  images 
to  compute  a  local  impulse  response  function  (LIR).  The  LIR  is  a  generalization  of  the  point-spread  function  applicable 
when  the  image  acquisition  and  reconstruction  processes  are  not  shift-invariant,  as  is  the  case  here.  We  then  computed  the 
two-dimensional  Fourier  transforms  of  the  LIRs  to  obtain  what  Barrett  and  Myers  have  called  a  local  modulation  transfer 
function  (LMTF)  [12].  A  standard  modulation  transfer  function  (MTF)  is  the  absolute  value  of  the  Fourier  transform  of 
the  system’s  point  spread  function  and  predicts  the  ratio  of  output  modulation  to  input  modulation  as  a  function  of  spatial 
frequency.  Localized  modulation  transfer  function  (LMTF)  is  the  generalization  of  MTF  to  linear,  shift-variant  systems. 
Higher  and  broader  LMTF  indicates  better  resolution  properties. 

Random  Gaussian  noise  with  mean  0  and  a  standard  deviation  of  1 .0  was  used  for  noisy  pressure  signals  for  the  noise 
studies.  Noisy  images  were  constructed  on  a  64  x  64  grid  using  a  zoom  factor  of  10.  The  noise  studies  were  performed  for 
500  realizations.  LNPS  is  a  generalization  of  the  noise  power  spectra  (NFS)  concept  that  can  be  used  for  linear  systems 
without  the  assumption  of  shift  invariance  which  does  not  hold  for  finite  transducer  aperture  systems  in  OAT.  LNPS  was 
computed  by  first  generating  a  set  of  500  realizations  of  reconstructed  images  for  each  algorithm  corresponding  to  pure 
Gaussian  noise  pressure.  For  each  set  of  these  reconstructed  images,  the  mean  image  was  computed.  The  mean  image  was 


then  subtracted  from  the  other  500  images.  LNPS  was  then  obtained  by  taking  the  average  squared  modulus  of  the  Fourier 
transform  (FT)  of  the  subtracted  images: 


LNPS 


\FT{noisyImage{i)  —  meanlmage)\^ , 


where  N  is  the  number  of  realizations. 


LNEQ  is  defined  as  the  ratio  of  the  square  of  the  LMTF  to  the  LNPS: 


LNEQ 


{LMTFf 

LNPS 


LNEQ  is  a  kind  of  frequency-dependent  signal-to-noise  ratio  generalized  to  linear,  shift-variant  systems.  Higher  LNEQ 
implies  higher  signal  detectability  performance  for  the  so-called  ideal  observer  in  the  task  when  both  signal  and  background 
are  known  exactly  [12]. 


3.  RESULTS 

3.1.  Phantom  Images 

In  all  the  phantom  images  shown  in  this  section,  the  transducer  axis  is  along  the  bottom  of  the  images.  Ligure  1  shows 
non-zoomed-in  images  produced  by  the  3  algorithms  for  a  circular  phantom  of  radius  1  mm  placed  at  a  distance  of  2  mm 
from  the  transducer  axis.  The  image  reconstructed  via  the  Lourier  based  algorithm  has  sharp  edges  but  it  is  non-uniform 
and  has  lower  contrast.  The  synthetic  aperture  image  is  quite  blurred.  The  image  reconstructed  via  Norton’s  approach 
is  quite  sharp  and  fairly  uniform.  Note  that  we  observe  two  images  in  the  Lourier-based  algorithm  due  to  the  implicit 
symmetry  assumption  in  the  reconstruction.  Ligure  2  shows  non-zoomed-in  images  of  a  line  of  small  rectangles  of  size  0.5 
mm  X  0.3  mm  placed  at  a  distance  of  2.0  mm  from  the  transducer  axis.  The  circular  arc  artifacts  are  more  visible  in  the 
synthetic  aperture  algorithm  than  the  Norton-based  algorithm  due  to  the  additional  filtration  step  that  is  performed  in  the 
Norton-based  algorithm.  In  addition,  the  rectangles  themselves  are  much  sharper  and  more  filled-in  in  the  Norton-based 
algorithm  compared  to  the  other  two  algorithms. 

3.2.  Spatial  Resolution 

The  images  of  zoomed-in  point  sources  are  shown  in  figure  3,  where  the  transducer  axis  is  along  the  bottom  of  the  images. 
The  LIR  plot  for  the  3  algorithms  is  shown  in  figure  4.  These  show  that  Lourier  based  algorithm  shows  the  best  resolution 
perpendicular  to  the  transducer  array  while  Norton’s  algorithm  shows  the  best  resolution  parallel  to  the  transducer  array. 
The  main  difference  between  Norton’s  algorithm  and  the  SA  algorithm  is  filtering.  This  results  in  a  much  narrower  LIR 
for  Norton’s  algorithm  than  for  the  SA  algorithm.  In  general,  the  lateral  resolution  for  Norton-based  and  SA  algorithms 
is  much  better  than  the  depth  resolution  (perpendicular  to  the  transducer  axis).  The  full  width  at  half-maxima  (LWHM) 
results  are  shown  in  table  1 . 

Ligure  5  shows  the  LMTL  images  for  the  three  algorithms.  The  LMTL  images  exhibit  an  asymmetry  due  to  the 
finite  transducer  length.  The  reciprocal  relationship  between  LIR  and  LMTL  is  exhibited  in  these  images.  LMTL  for  SA 
algorithm  is  the  narrowest  since  LIR  for  SA  is  the  broadest.  Ligure  6  shows  the  LMTL  plots.  LMTL  for  Lourier  based 
algorithm  is  the  best  in  the  direction  perpendicular  to  the  transducer  axis,  especially  for  smaller  frequencies.  Norton’s 
Algorithm  produces  the  best  LMTL  profile  in  the  lateral  direction  which  is  expected  since  it  had  the  smallest  lateral 
resolution. 


3.3.  Noise  texture 


The  noise  texture  images  are  simply  images  of  a  uniform  background  reconstructed  from  pressure  signals  to  which  zero 
mean  Gaussian  noise  has  been  added.  Figure  7  shows  the  noise  texture  in  the  unzoomed  images  reconstructed  via  the 
three  algorithms.  The  noise  in  the  Fourier  and  Norton-based  algorithms  seems  uniformly  speckled  while  the  smeared  out 
noise  texture  in  SA  algorithm  seems  to  exhibit  some  long  range  correlations.  Such  ’blobbiness’  in  the  noise  can  impede 
detectability  of  signals  of  size  comparable  to  the  blob  size.  LNPS  describes  the  frequency  content  of  the  background 
texture  in  the  region  surrounding  the  location  of  the  signal  in  the  image  [12].  Figure  8  shows  the  zoomed-in  (250%) 
images  of  LNPS  for  the  three  algorithms.  LNPS  images  for  Norton-based  and  synthetic  aperture  algorithms  are  pretty 
symmetric.  But  it  is  not  so  for  Fourier-based  algorithm.  Note  that  the  input  to  the  Fourier-based  algorithm  was  Gaussian 
noise  pressure  while  the  input  to  the  other  two  algorithms  was  time-integrated  Gaussian  noise  pressure,  which  introduces 
noise  correlations  that  can  affect  the  form  of  the  LNPS.  Figure  9  shows  the  LNPS  plots  in  the  two  directions.  The  plot  for 
SA  algorithm  was  omitted  because  it  was  several  orders  of  magnitude  higher  than  the  other  two  algorithms  and  could  not 
be  fit  into  the  same  plot.  The  shapes  for  LNPS  are  very  different  for  Norton-based  and  Fourier-based  algorithms. 

3.4.  Signal  Detectability/LNEQ 

Figure  10  shows  the  zoomed-in  (250%)  images  of  LNEQ  for  the  three  algorithms.  Figure  11  shows  the  LNEQ  plots. 
In  general,  LNEQ  images  are  somewhat  difficult  to  interpret  visually,  as  they  represent  a  kind  of  generalized  frequency- 
domain  detectability  transfer  function,  but  bright  areas  correspond  to  frequency  components  that  are  more  likely  to  be 
detected  against  a  background  of  correlated  noise,  as  captured  by  the  LNPS.  To  calculate  the  so-called  ideal  observer 
signal  to  noise  ratio  for  a  small  low  contrast  signal,  one  would  calculate  the  squared  magnitude  of  the  Eourier  transform 
of  the  signal,  multiply  it  by  the  LNEQ  shown  in  the  images,  and  integrate  over  all  frequencies.  The  plots  in  Eig.  11 
give  a  better  sense  of  the  relative  magnitude  of  the  LNEQ  for  the  different  algorithms.  Higher  values  are  better  and  the 
Norton-based  algorithm  produces  the  highest  LNEQ  in  both  directions  followed  by  the  Eourier-based  algorithm.  This 
indicates  superior  ideal  observer  signal  detectability  in  images  reconstructed  by  the  use  of  the  Norton-based  algorithm.  Eor 
all  algorithms  the  LNEQ  becomes  small  at  5  mm  '  in  the  depth  direction  and  7  mm  '  in  the  lateral  direction. 

4.  DISCUSSION 

The  2D  geometry  that  we  consider  in  this  paper  reduces  the  inherently  3D  optoacoustic  image  reconstruction  to  2D  due  to 
transducer’s  2D  directivity.  This  results  in  the  measured  2D  optoacoustic  signal  being  related  to  the  time  derivative  of  the 
integral  of  the  absorbed  optical  energy  over  circular  arcs  centered  at  the  transducer  as  given  by  eqn.  3.  The  three  algorithms 
considered  here  for  evaluation  have  some  intrinsic  differences.  The  SA  algorithm  is  approximate  to  start  with  and  involves 
only  the  delay  and  sum  technique  based  on  backprojection  over  circular  arcs.  The  Norton-based  algorithm  improves  on 
SA  algorithm  by  providing  a  filtration  step  that  filters  the  time-integrated  pressure  data  before  using  the  backprojection 
over  circular  arcs  technique.  The  Eourier-based  algorithm  is  not  exact  in  the  2D  geometry  that  we  are  considering  since 
it  assumes  that  the  optoacoustic  source  is  constant  in  the  third  dimension.  These  intrinsic  differences  in  the  algorithms 
explain  the  differences  in  sharpness  and  quality  of  the  reconstructed  images. 

Norton-based  and  SA  algorithms  use  time-integrated  pressure  as  input  signal  while  the  Eourier-based  algorithm  uses 
pressure  signal  as  input.  Integration  of  noisy  data  introduces  noise  correlations  that  can  affect  the  LNPS.  This  was  reflected 
in  the  blobbiness  of  the  noise  texture  images  of  the  SA  algorithm.  The  additional  filtration  step  in  the  Norton-based 
algorithm  aids  in  the  removal  of  such  blobbiness  and  gives  a  much  better  noise  texture.  The  Eourier-based  algorithm  does 
not  show  such  blobbiness  in  noise  texture. 


5.  CONCLUSIONS 

We  explored  three  different  ways  in  which  a  2D  image  can  be  reconstructed  in  OAT.  We  implemented  and  evaluated  three 
algorithms  for  2D  image  reconstruction  in  OAT-  Eourier-based,  Norton-based  and  synthetic  aperture  algorithms.  We  found 
that  the  2D  Eourier-based  algorithm  offers  better  resolution  and  LMTE  in  the  depth  direction  while  Norton’s  algorithm 
offers  the  best  lateral  resolution.  However,  we  found  that  in  reconstructions  of  phantoms  the  images  produced  by  Norton- 
based  algorithm  looked  the  sharpest  and  more  uniform.  The  LNEQ  data  suggests  that  Norton-based  algorithm  has  the  best 
signal  detectability. 
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6.  APPENDIX 

A.  2D  relation  between  pressure  and  absorbed  optical  energy  function 

Starting  with  equation  (2), 


p{r,t)  =  77 


d^r'A(r') 


d  S{t  - 


dt  47r|r 


C 


=  /  Aflds' 

47r  dt  J  ct 

I  Ar— ci| 


(18) 


where  Ar  =  |r  —  r'|  ,  dS'  is  the  differential  surface  area  and  the  integral  is  over  a  spherical  surface  centered  on  r  and  of 
radius  Ar.  This  can  be  written  as: 


I  Ar|=ct 

where  dO'  =  sin  9'd9'd(j)'  is  the  solid  angle.  This  equation  in  2D  (x  —  z  plane)  reduces  to: 

P(a^,  2,  f)  =  ^  J  {ct)A{p')d<j)' , 

\p-p'\=ct 

where  p  =  \/x‘^  +  ,  is  the  polar  radial  coordinate. 

Using  the  integral  property  of  delta  functions,  this  can  be  written  in  polar  coordinates  as: 

p{x,z,t)  =  J  J  \p- p'\d{\p- p'\)S{ct-\p- p'\)A{p')d(l)\ 

\p-p>\=ct 


(19) 


(20) 


This  can  be  written  equivalently  in  2D  Cartesian  coordinates  as: 

p{x,z,t)  j  j  dx'dz' A{x' ,  z')d  —  \J {x  —  x'Y  +  . 

B.  Derivation  of  Norton-based  algorithm 

In  this  section,  we  shall  derive  Eq.  (12)  following  the  method  detailed  in  Norton’s  paper  [11].  Define  the  Eourier  transform 
with  respect  to  r  as: 


f{x,v) 


f{x,  r)  exp{i2T:vr)dr. 


Taking  the  Eourier  transform  of  equation  (10)  on  both  sides  with  respect  to  p  =  we  get 


(21) 


g'{x^  v)  =  A' [x,  v)  *  exp(727ra:^77). 


(22) 


On  convolving  both  sides  with  exp(— we  get 


g'{x,  v)  *  exp(— i27ra;^^)  =  A' {x,  v)  *  exp(i27ra;^i^)  *  exp(— z27rx^i^), 
where  the  convolution  is  with  respect  to  x.  Using  the  identity 


exp(z27ra;^zz)  *  exp(— z27ra;^zz)  = 


5{2ux) 

6{x)  6{v) 


2\v\  2|x 


equation  (23)  becomes 


g'{x,  v)  *  exp(— z27ra;^zz)  =  A'{x,  v)  *  (  — ^  + 


()(a:)  5{v) 


2W\  2x 


Solving  for  A'{x,  v)  gives, 


(23) 


(24) 


A'{x,  v) 


2\v\g'{x,  v)  *  exp(— z27ra;^j^)  —  \v\5{v) 


Using  the  identity  \v\5{v)  =  0  to  eliminate  the  second  term  on  the  right  and  taking  the  inverse  Fourier  transform 
(FT“^)with  respect  to  v  on  both  sides  one  finds 


A'{x,  p)  =  g'{x,  p)  *  *FT  ^{2\v\ exp(— z27ra;^j^)}p. 
Comparing  Eq.  (25)  with  Eq.(l  1),  we  see  that 


(25) 


R{x,p)  =  FT  ^{2|zz|  exp(— z27ra;^zz)}p 

=  2FT-i{k|}p+.2  (26) 

If  the  data  S'x,  v)  is  bandlimited  by  a  cut-off  frequency  then  the  above  convolution  relation  is  unchanged  if  we  impose 
the  same  bandlimit  on  R{x,  p).  So  one  gets: 


R{x,  p)  =  2FT  <  \iy\vect  — 


2Vr 


One  can  write  out  the  convolution  in  equation  (25)  as 


A'{x,p)  =  vl  /  I  g'{x',p)Rx[vAC- P+i.x-x'Y)]dpdx', 


where  Ri{u)  =  4sinc(2zz)  —  2sinc^(zz). 

Substituting  for  A' {x,  p),  S'{x,  p)  and  for  (  and  p,  one  obtains 


A{x,  z)  =  2zv'l  r)_Ri[zzc(^^  —  +  {x  —  x'Y)]drdx' , 


which  is  the  desired  result. 


(27) 


(28) 
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FWHM 

Fourier-based 

Norton-based 

Synthetic  Aperture 

Depth  (mm) 

0.154 

0.200 

0.471 

Lateral  (mm) 

0.161 

0.151 

0.189 

Table  1.  FWHM  values  for  a  point  source  of  size  0.1  mm  with  pixel  width=0.1  mm 


(a)  Fourier  -based  (b)  Norton-based  (c)  Synthetic  Aperture 

Figure  1.  Circular  phantom  images 


(a)  Fourier-based  (b)  Norton-based  (c)  Synthetic  Aperture 

Figure  2.  Images  of  a  line  of  rectangles 


(a)  Fourier-based  (b)  Norton-based  (c)  Synthetic  Aperture 


Figure  3.  Zoomed-in  point  source  images 


(a)  Fourier-based 


(b)  Norton-based 


(c)  SA 


Figure  5.  Zoomed-in  LMTF  Images 


(a)  Perpendicular  to  the  transducer  axis  (b)  Parallel  to  the  transducer  axis 


Figure  6.  LMTF  plots 


(a)  Fourier-based 


(b)  Norton-based 


Figure  8.  Zoomed-in  LNPS  Images 


(a)  Perpendicular  to  the  transducer  axis 


(b)  Parallel  to  the  transducer  axis 


Figure  9.  LNPS  Plots 


(a)  Fourier-based  (b)  Norton-based  (c)  SA 


Figure  10.  LNEQ  Images 
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ABSTRACT 

Previous  research  correcting  for  variable  speed  of  sound  in  photoacoustic  tomography  (PAT)  has  used  a  generalized  radon 
transform  (GRT)  model .  In  this  model,  the  pressure  is  related  to  the  optical  absorption,  in  an  acoustically  inhomogeneous 
medium,  through  integration  over  non-spherical  isochronous  surfaces.  This  model  assumes  that  the  path  taken  by  acoustic 
rays  is  linear  and  neglects  amplitude  perturbations  to  the  measured  pressure.  We  have  derived  a  higher-order  geometrical 
acoustics  (GA)  expression,  which  takes  into  account  the  first-order  effect  in  the  amplitude  of  the  measured  signal  and 
higher-order  perturbation  to  the  travel  times.  The  higher-order  perturbation  to  travel  time  incorporates  the  effect  of  ray 
bending.  Incorrect  travel  times  can  lead  to  image  distortion  and  blurring.  These  corrections  are  expected  to  impact  image 
quality  and  quantitative  PAT.  We  have  previously  shown  that  travel-time  corrections  in  2D  suggest  that  perceivable  differ¬ 
ences  in  the  isochronous  surfaces  can  be  seen  when  the  second-order  travel-time  perturbations  are  taken  into  account  with 
a  10%  speed  of  sound  variation.  In  this  work,  we  develop  iterative  image  reconstruction  algorithms  that  incorporate  this 
higher-order  GA  approximation  assuming  that  the  speed  of  sound  map  is  known.  We  evaluate  the  effect  of  higher-order 
GA  approximation  on  image  quality  and  accuracy. 

Keywords:  Optoacoustic  tomography,  photoacoustic  tomography,  thermoacoustic  tomography,  image  reconstruction,  vari¬ 
able  speed  of  sound,  inhomogeneous,  travel  times,  geometrical  acoustics. 


1.  INTRODUCTION 

Optoacoustic  imaging  is  a  hybrid  imaging  technique  that  has  attracted  a  lot  of  attention  in  recent  years.  It  is  based  on  the 
photoacoustic/optoacoustic  effect  which  refers  to  acoustic  wave  generation  upon  absorption  of  pulsed  optical  energy  by  a 
medium.  A  slight  rise  in  temperature  of  the  medium  due  to  the  absorption  of  the  incident  electromagnetic  wave  results  in 
thermoelastic  expansion.  This  thermoelastic  expansion  and  then  contraction  due  to  the  pulsed  electromagnetic  waves  leads 
to  the  generation  of  acoustic  waves.  Under  the  constraints  of  thermal  and  stress  confinement,  this  thermal  expansion  leads 
to  a  rise  in  pressure,  p{r,t),  that  satisfies  the  three-dimensional  inhomogeneous  wave  equation  [1]: 


df^ 


OpWt 


(1) 


where  the  heating  function,  is  the  thermal  energy  deposited  by  the  electromagnetic  radiation  per  unit  time  per  unit 

volume,  [3  is  the  isobaric  volume  expansion  coefficient  and  Cp  is  the  specific  heat  of  the  medium.  The  heating  function  can 
be  expressed  as  the  product  of  a  spatially  varying  optical  absorption  function  of  the  medium,  A(r),  and  a  time  dependent 
optical  illumination  function,  / (f).  Using  standard  Green’s  function  techniques,  the  measure  pressure  signal  can  be  related 
to  the  optical  absorption  function,  assuming  delta  pulse  illumination,  as: 


p{r,t)  =  T] 


d^r'A{r') 


d  S{t- 


dt  47r|r 
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(2) 


where  77  =  ^ .  Eqn.  (2)  states  that  the  time  integral  of  acoustic  pressure  at  a  point  r  and  time  t  is  given  by  the  integral 
of  the  optical  absorption  function  over  a  spherical  surface  of  radius  |r  —  r'|  =  ct  centered  at  r.  A  simple  but  inexact 
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way  to  reconstruct  A(r)  is  to  spatially  resolve  the  optoacoustic  waves  by  using  the  speed  of  sound  and  to  backproject 
time-integrated  pressure  signals  over  hemispheres. 

This  backprojection  assumes  that  the  speed  of  sound  in  the  medium  is  constant  and  known.  However,  both  these 
assumptions  are  not  valid  in  tissues.  Speed  of  sound  in  tissues  can  vary  from  1350  m/s  (in  fat)  to  about  1700  m/s  (for 
skin)  for  ultrasound  waves  in  the  1-10  MHz  range  [2].  Using  a  constant  speed  of  sound  in  image  reconstruction  can  result 
in  image  distortion  and  poor  image  quality.  There  has  been  some  work  done  previously  to  address  the  effect  of  speed  of 
sound  variations  in  OAT.  Xu  and  Wang  [3]  [4]  looked  at  variable  speed  of  sound  reconstruction  in  breast  thermoacoustic 
tomography  (TAT).  They  concluded  that  there  is  minor  amplitude  distortion  in  breast  TAT  and  that  variation  in  travel  time 
is  a  first-order  perturbation  in  a  weakly  acoustically  heterogeneous  medium.  Zhang  and  Anastasio  [5]  derived  a  heuristic 
method  for  reconstructing  acoustic  speed  and  optical  absorption  distributions  in  a  step-wise  manner.  Their  work  was 
also  based  on  the  first-order  travel  time  effect.  Kolkman  et  al.  [6]  devised  a  method  to  determine  the  speed  of  sound  in 
tissue  using  two  concentric  rings  based  photoacoustic  sensor  based  on  first-order  travel  time  effect.  Manohar  et  al.  [7] 
also  devised  a  new  and  improved  method  to  determine  speed  of  sound  using  a  photoacoustic  imager.  Agravonsky  and 
Kuchment  [8]  have  recently  derived  an  analytic  reconstruction  formula  for  OAT  for  arbitrary  detector  geometry,  such 
that  the  point  detectors  are  placed  on  a  closed  surface  and  for  variable  speed  of  sound.  Their  analytic  formula  leads  to 
reconstruction  in  terms  of  eigenfunction  expansion.  Zhen  et  al.  [9]  devised  an  iterative  reconstruction  algorithm,  using 
finite-element  method,  that  incorporates  both  attenuation  and  variable  speed  of  sound. 

Previous  work  on  the  speed  of  sound  variation  that  is  based  on  the  generalized  radon  transform  (GRT)  model  has  only 
looked  at  the  first-order  effect  of  variable  speed  of  sound  [4]  [5]  [6].  In  the  GRT  model  the  pressure  is  given  by: 


p{r,t)  =  r] 


d\'A{r') 


d  S{t  —  T(r,  r')) 
dt  47r|r  — r'l  ’ 


(3) 


where  T(r,r')  is  the  travel  time  for  the  sound  wave  to  travel  between  points  r  and  r'.  It  was  assumed  that  sound  rays 
continue  to  travel  in  straight  line  in  the  presence  of  acoustic  heterogeneity  and  the  travel-time  between  points  r  and  r'  is 
given  as  a  line  integral  over  the  slowness  map: 


T(r,r')  =  J  ds/c{s),  (4) 

ro 

where  Tq  is  the  straight  line  joining  the  points  r  and  r',  ds  is  the  differential  length  element  along  that  line  and  c(s)  is  the 
variable  speed  of  sound.  The  effect  of  variable  speed  of  sound  on  signal  amplitude  was  neglected.  Previous  work  based  on 
the  GRT  model  does  not  incorporate  the  effect  of  ray  bending  which  may  be  significant  when  the  speed  of  sound  varies  by 
10%  or  more. 

In  this  work,  we  investigate  a  higher-order  geometrical  acoustic  approximation  that  incorporate  first-order  effect  on  the 
signal  amplitude  and  higher-order  effect  on  the  travel  times.  We  use  this  higher  order  approximation  to  construct  the  system 
matrix.  We  then  iteratively  reconstruct  the  images  using  third-order,  second-order,  first-order  and  zeroth-order  travel  time 
effects.  We  show  that  the  higher-order  GA  approximation  offers  much  better  image  quality  and  accuracy  when  the  speed 
of  sound  varies  by  10%  or  more. 


2.  METHODS 

We  treat  the  variation  in  the  speed  of  sound  as  a  perturbation  to  the  background  sound  speed,  cq. 

c(r)  =  Co  +  eci(r), 

where  e  characterizes  the  magnitude  of  the  perturbation. 


(5) 


2.1.  Derivation  of  Geometrical  Acoustics  Approximation 

The  GA  approximation  ignores  diffraction  effects.  It  is  valid  in  the  short  wavelength  regime  when  the  size  of  the  inhomo¬ 
geneity  is  much  greater  than  the  wavelength.  In  a  scattering  medium,  this  approximation  is  valid  when  the  speed  of  sound 
does  not  change  significantly  over  one  wavelength.  The  Green’s  function  in  the  acoustically  inhomogeneous  medium  can 
be  written  as: 


G(r,r',w)  =  5o(r,r')exp(iwr(r,r')),  (6) 

where  t  is  known  as  the  eikonal  function.  The  GA  approximation  gives  us  the  eikonal  equation  in  the  limit  A  — >  0: 

[VrT(r,r')]  =  c"2(r),  (7) 

The  eikonal  equation  implies: 

r 

=  J  ds/c{s),  (8) 

r' 

where  s  denotes  the  arc  length  along  the  path  between  the  points  r  and  r'  as  shown  in  figure  1.  Note  that  the  curve 
connecting  the  points  r  and  r'  need  not  be  a  straight  line.  The  ray  trajectory  is  determined  by  the  path  that  minimizes  the 
acoustic  path  length  (Fermat’s  principle)  or  equivalently,  that  minimizes  the  travel-time.  For  an  acoustically  homogeneous 
medium,  this  is  a  straight  line.  Previous  work  by  Snieder  et  al.  [10]  concludes  that  to  first-order  in  perturbation,  this 
trajectory  can  chosen  to  be  the  path  along  the  reference  ray  that  satisfies  the  eikonal  equation  (assuming  the  variation  in 
speed  of  sound  is  slowly  varying). 


Figure  1.  The  path  travelled  by  sound  wave  between  points  r  and  r’. 


Using  Eqn.  (6),  the  pressure  is  given  by: 


p{v,oj)  ~ —iujrj  J  d^r'A{r')go{r,r')  exp(iwT(r,  r')), 

Taking  Fourier  transform  with  respect  to  uj  we  obtain  the  generalized  radon  transform  (GRT): 

p(r,f)~77  J  d^r'A{A)go{r,r')^d{t-T{r,r')). 

To  zeroth-order,  the  amplitude  is  given  by: 

Using  this  expression  for  5o(r)  rO’  we’ll  obtain  the  following  expression  for  pressure: 


p(r,  0  ,  I  {t  -  T(r,  r'))  ■ 


(9) 


(10) 


(11) 


(12) 


This  expression  along  with  the  travel  time  given  by  Eqn.  (4)  is  what  has  been  used  in  the  previous  work  addressing 
speed  of  sound  variation  based  on  the  GRT  model. 


2.2.  Higher-order  Geometrical  Acoustics  Approximation 

We  can  improve  the  first-order  GA  model  by  incorporating  higher  order  effects  on  the  amplitude  and  travel  times. 


•  Improvements  to  amplitude  -  The  first-order  correction  to  the  amplitude  of  the  Green’s  function  is  given  by: 


ffo(r,r') 


1  /  c(r) 

47r|r  —  r'l  y  c(r') 


(13) 


•  Improvements  to  the  eikonal  -  We  can  incorporate  the  effect  of  ray  bending  by  considering  the  higher  order  pertur¬ 
bations  in  the  eikonal.  The  assumption  for  first-order  GA  is  that  the  speed  of  sound  is  slowly  varying  so  that  the  time 
of  travel  can  be  obtained  using  linear  rays.  However,  if  this  assumption  is  not  true,  it  can  lead  to  higher-order  pertur¬ 
bations  in  travel  times.  Higher-order  travel  time  perturbations  contribute  to  ray  bending.  Following  the  methodology 
of  Snieder  et  al.  [10],  perturbations  in  travel-time  can  be  written  as  : 


T(r,r')  =  To(r,r')  -h  eTi(r,r')  -h  e^r2(r,r')  -h  .....  (14) 

where  e  is  a  small  parameter. 

Let  To  denote  the  reference  ray  associated  with  the  reference  eikonal  ro(r,r')  =  ds/cQ.  Let  the  reference  ray  be 
parametrized  by  the  variable  s  so  that  Tq  =  sto,  where  to  is  the  unit  vector  along  the  reference  ray.  The  first-  and  second- 
order  perturbations  to  travel  time  are  given  by  : 


Ti(r,r')  =  -  [  (15) 

Jro  Cq 

and 

T2(r,r')=^  ^  -  |Vri(r,r')n^  ds  (16) 

The  higher-order  perturbations  to  travel  times  can  be  easily  derived  as  described  by  Snieder  et  al.  [10].  Thus,  one 
can  calculate  the  perturbed  time  of  flight  using  these  equations  if  one  knows  the  reference  ray.  Note  that  this  perturbation 
theory  approach  only  works  if  the  non-linear  perturbations  are  small  and  there  is  no  multipathing.  Thus,  we  can  obtain 
a  higher-order  estimate  of  the  pressure  than  that  given  by  Eqn.  (12)  by  using  Eqn.  (13)  for  amplitude  and  keeping  up  to 
higher-order  terms  in  Eqn.  (14)  and  using  these  in  Eqn.  (10). 

2.3.  Details  of  the  simulation 
2.3.1.  Travel  time  calculations 

Travel  times  were  calculated  using  up  to  fourth-order  correction  for  different  pixels  for  a  specific  transducer  position.  This 
was  done  in  2D  for  a  50x50  grid  for  a  specific  speed  of  sound  map  shown  in  figure  2.  The  blurred  elliptical  region  (2  mm 
by  1 .4  mm)  in  this  map  has  the  variable  speed  of  sound  with  respect  to  the  background.  The  background  speed  of  sound 
was  set  to  1500  m/s.  Curves  of  constant  times  were  reconstructed  using  up  to  second-order  travel  times.  This  was  done  for 
a  5%,  10%  and  15%  speed  of  sound  variation.  The  pixel  size  was  set  to  0.01  cm.  The  time  sampling  interval  was  set  to 
33.33  ns.  The  transducers  were  assumed  to  lie  on  a  line  to  the  left  of  the  phantom  0.01  cm  apart. 


Figure  2.  Speed  of  sound  map 


2.3.2.  Image  Reconstruction 

The  system  matrices  were  constructed  in  2D  incorporating  zeroth-order,  first-order,  second-order  and  fourth-order  travel 
time  effects.  The  images  were  then  reconstructed  iteratively  using  least-squares  method.  The  conjugate  gradient  method, 
was  used  to  choose  the  step  direction  for  each  iteration.  The  fourth-order  travel  time  system  matrix  was  used  for  the 
forward  model.  The  images  were  then  reconstructed  iteratively  using  third-order,  second-order,  first-order  and  zeroth- 
order  matrices.  The  image  reconstruction  was  done  on  a  50x50  grid  with  parameters  as  specified  above. 

3.  RESULTS 

3.1.  Isochrons  and  relative  travel  times 

Figure  3  depicts  the  change  to  relative  travel  times  for  pixels  located  along  the  line  (0,  1 .0  mm)  with  the  transducer  placed 
at  (0,  2.3  mm)  in  pixel  coordinates.  One  observes  that  the  higher-order  travel  time  corrections  become  perceivable  when 
the  speed  of  sound  varies  by  10%  or  more.  One  can  also  see  that  the  travel  time  perturbations  seem  to  converge  when  one 
goes  up  to  the  third-order  term.  This  suggests  that  one  should  use  up  to  third-order  travel  time  corrections  in  the  forward 
model  to  most  accurately  represent  the  signal. 


Relative  travel  limes<5%  variation) 


Relative  travel  times<10%  variation) 


Relative  travel  limesll§%  variation) 


Pixel# 


Pixel# 


(a)  5  percent  variation 


(b)  10  percent  variation 


(c)  15  percent  variation 


Figure  3.  Relative  travel  times  for  transducer  at  (0,  2.3  mm)  for  phantom  pixels  along  the  line  (0,1.0  mm)  for  a  specific  speed  of  sound 
map 


Figures  4  and  5  depict  the  isochrons  for  the  speed  of  sound  map  shown  in  figure  2  with  the  transducer  placed  at  (0,  2.0 
mm)  and  for  a  travel  time  of  1.665  micro-seconds. 


(a)  Zeroth-order  (b)  First-order  (c)  Second-order 
Figure  4.  Isochrons  for  transducer  at  (0,  2.0  mm)  for  time=166.65  -seconds  with  a  10%  speed  of  sound  variation 


(a)  Zeroth-order  (b)  First-order  (c)  Second-order 

Figure  5.  Isochrons  for  transducer  at  (0,  2.0  mm)  for  time=166.65  nano-seconds  with  a  15%  speed  of  sound  variation 


From  these  figures,  it  can  be  concluded  that  a  speed  of  variation  of  10%  or  more  leads  to  some  differences  between  the 
isochrons  incorporating  second-order  and  first-order  travel  time  effects. 

3.2.  Iteratively  reconstructed  images 

Figures  6  and  7  show  the  iteratively  reconstructed  images  incorporating  up  to  third-,  second-,  first-  and  zeroth-order  travel 
time  effects.  For  a  5%  speed  of  sound  variation,  the  images  reconstructed  using  third-  and  second-order  corrections  are 
pretty  close  to  the  true  phantom.  The  first-order  correction  offers  a  better  image  than  the  zeroth-order  correction,  however, 
it  is  is  not  as  good  as  the  images  reconstructed  using  higher-order  corrections.  For  a  10%  speed  of  sound  variation,  the 
image  reconstructed  using  third-order  correction  is  much  closer  to  the  true  phantom  than  those  reconstructed  using  the 
lower-order  corrections. 


(a)  True  (b)  Third-order  (c)  Second-order  (d)  First-order  (e)  Zeroth-order 

Figure  6.  Iteratively  reconstructed  images  using  fourth-order  travel  time  correction  in  the  forward  model  with  a  5%  speed  of  sound 
variation,  Image  were  reconstructed  using  various-order  travel  time  corrections 


(a)  True  (b)  Third-order  (c)  Second-order  (d)  First-order  (e)  Zeroth-order 

Figure  7.  Iteratively  reconstructed  images  using  fourth-order  travel  time  correction  in  the  forward  model  with  a  10%  speed  of  sound 
variation,  Image  were  reconstructed  using  various-order  travel  time  corrections 


Figure  8  shows  a  line  profile  through  the  phantom  at  (0,  2.3  mm).  One  can  see  that  the  difference  between  third-order 
and  second-order  is  not  so  much  for  a  5%  speed  of  sound  variation.  However,  the  difference  is  more  significant  for  a  10% 
speed  of  sound  variation  especially  for  objects  lying  behind  the  inhomogeneity.  The  third-order  correction  offers  the  best 
match  to  the  true  phantom  in  this  case. 


(a)  5  percent  variation  (b)  10  percent  variation 

Figure  8.  Line  profile  through  the  phantom  at  (0,2.3  mm) 


4.  CONCLUSIONS 

We  derived  a  higher-order  geometrical  acoustics  approximation  to  the  GRT  model  in  PAT.  This  incorporates  the  first-order 
correction  to  the  pressure  amplitude  and  higher-order  correction  to  the  travel  times.  We  found  that  some  differences  can 
be  seen  in  the  isochronous  curves  between  second-  and  first-order  corrections  when  the  speed  of  sound  varies  by  10% 
or  more.  These  differences  in  travel  times  translated  into  the  reconstructed  images  as  well.  Images  that  were  iteratively 
reconstructed  using  the  third-order  travel  time  corrections  were  more  accurate  than  those  constructed  using  the  lower-order 
corrections.  We  will  conduct  further  studies  using  actual  pressure  data  in  an  inhomogeneous  medium  in  3D  to  quantify  the 
effect  of  higher-order  GA  approximation. 
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ABSTRACT 

Attenuation  effects  can  be  significant  in  photoacoustic  tomography  (PAT)  since  the  measured  pressure  signals  are  broad¬ 
band  and  ignoring  them  may  lead  to  image  artifacts  and  blurring.  Previous  work  by  our  group  had  derived  a  method  for 
modeling  the  attenuation  effect  and  correcting  for  it  in  the  image  reconstruction.  This  was  done  by  relating  the  ideal, 
unattenuated  pressure  signals  to  the  attenuated  pressure  signals  via  an  integral  operator.  In  this  work,  we  explore  singular- 
value  decomposition  (SVD)  of  a  previously  derived  3D  integral  equation  that  relates  the  Fourier  transform  of  the  measured 
pressure  with  respect  to  time  and  two  spatial  components  to  the  2D  spatial  Fourier  transform  of  the  optical  absorption 
function.  We  find  that  the  smallest  singular  values  correspond  to  wavelet-like  eigenvectors  in  which  most  of  the  energy  is 
concentrated  at  times  corresponding  to  greater  depths  in  tissue.  This  allows  us  characterize  the  ill  posedness  of  recovering 
absorption  information  from  depth  in  an  attenuating  medium.  This  integral  equation  can  be  inverted  using  standard  SVD 
methods  and  the  optical  absorption  function  can  be  recovered.  We  will  conduct  simulations  and  derive  algorithm  for  image 
reconstruction  using  SVD  of  this  integral  operator. 

Keywords:  Optoacoustic  tomography,  photoacoustic  tomography,  thermoacoustic  tomography,  image  reconstruction,  at¬ 
tenuation,  singular-value  decomposition 


1,  INTRODUCTION 

Optoacoustic/photoacoustic  imaging  is  an  emerging  imaging  technique  [1,2]  based  on  the  photoacoustic/optoacoustic 
effect.  This  effect  refers  to  acoustic  wave  generation  upon  absorption  of  pulsed  optical  energy  by  a  medium.  The  majority 
of  the  image  reconstruction  algorithms  in  PAT  so  far  have  assumed  a  non-dispersive  acoustic  medium.  The  effect  of 
frequency-dependent  attenuation  on  acoustic  waves  can  be  significant  since  PAT  uses  broadband  detection.  Reconstructed 
images  may  exhibit  distortion  and  artifacts  if  these  effects  are  not  taken  into  account.  Previous  work  on  dispersive  acoustic 
media  done  by  our  group  [3]  focused  on  incorporating  the  frequency-dependent  attenuation  effects  into  the  PAT  model. 
This  was  done  by  using  a  ID  integral  operator  to  relate  the  ideal  pressure  to  the  attenuated  pressure.  We  have  also  previously 
derived  an  inversion  formula  for  the  optical  absorption  function  using  singular-value  decomposition  (SVD)  [4]  based  on 
an  approach  similar  to  that  of  Schotland  et  al.  [5]  [6],  in  optical  diffusion  tomography.  This  expression  directly  relates  the 
measured  attenuated  pressure  to  the  optical  absorption  function.  This  formula  is  applicable  in  a  planar  geometry  where 
the  array  of  transducers  lies  in  a  plane.  In  this  paper,  we  will  use  this  approach  for  image  reconstruction  in  an  attenuating 
medium.  We  will  also  conduct  noise  and  resolution  studies  on  the  images  reconstructed  using  this  approach. 

2.  METHODS 

In  an  attenuating  medium,  the  optoacoustic  wave  equation  includes  a  loss  term  [3]  and  is  given  by: 

where  L{t)  =  ^  J  dcu  ~ 

=  -^  +  iaiuj) ,  (2) 

c(w) 
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where  u;  is  the  temporal  frequency  and  k{uj)  is  the  complex  wave  number  [7]. 
For  tissues,  the  ultrasonic  attenuation  a  is  given  to  a  good  approximation  by: 


a{u})  =  ao|w|,  (3) 

where  cm‘'rad'^s  . 

Using  Kramers-Kronig  relations,  one  obtains  the  dependence  of  phase  velocity  c(u;)  on  frequency  due  to  frequency- 
dependent  attenuation  as: 


1  1  2  a; 

— -  = - aoln\  —  I, 

c(a;)  Co  TT  a;o 

where  a;o  is  the  reference  frequency  at  which  c(a;)  =  cq. 

So  for  an  attenuating  medium,  on  taking  the  temporal  Fourier  transform  of  Eqn.  (1)  we  get: 


(4) 


V^p(r,a;)  -f  k^{uj)p{v,i^)  =  - A{r) I (u) . 

Op 

This  equation  can  be  solved  using  standard  Green’s  function  methods  as  [8]: 


(5) 


p(r,a;)  =  —ir]U!l{u!) 


where  r]  =  ^  and  G(r  —  r') 


|r— r’  |  ^ 
47r|r— r'l 


d^r'A{r')G{r-r'), 


(6) 


2.1.  SVD  of  integral  operator  relating  attenuated  pressure  to  optical  absorption  function 

Using  the  angular  spectrum  expansion  of  Green’s  function  [9]  and  by  simplifying  and  discretizing  Eqn.  (6),  one  can  obtain 
an  integral  equation  [4] : 


-Pn(q) 


J  dzA{q,  z)iT„(q,  z). 


(7) 


This  equation  can  be  inverted  as  [6]: 


where 


and 


4(q,  ^)  =  51  z)M^l{q)Pr^{q), 

m,n 

Mmn{q)=  J  dzKm{q,z)K*{q,z). 


(8) 

(9) 

(10) 


The  two-dimensional  spatial  wave  vector  q  is  defined  as: 

q={k^,ky,0),  (11) 

and  the  temporal  frequency  oj  is  discretized  as  a;„,  F’n(q)  =  p(q,  w„)  and  =  fc(a;„)  . 

The  computation  of  the  inverse  of  this  matrix  M  is  the  key  step  in  the  procedure  to  recover  the  optical  absorption 
function.  Its  pseudoinverse  can  be  computed  by  performing  an  SVD  of  matrix  K.  Some  of  the  eigenvalues  of  this  matrix 
may  be  zero  or  very  small.  This  problem  can  be  avoided  by  using  regularization. 

Thus,  the  procedure  for  recovering  A(r)  is: 


1 .  Take  the  Fourier  transform  of  the  measured  pressure  (at  z  =  0)  with  respect  to  time  and  the  2D  spatial  components 
to  obtain  p{kx,ky,uj). 

2.  Discretize  p{kx,  ky,uj)  as  p(q,  w„). 

3.  Compute  the  matrix  Kn{q,  z)  using  Eqn.  (9)  . 

4.  Compute  the  pseudoinverse  of  matrix  Mmn{q)  by  performing  SVD  of  matrix  K. 

5.  Compute  A(q,  z)  using  Eqn.  (8). 

6.  Take  the  2D  inverse  Eourier  transform  of  A(q,  z)  to  obtain  A(r). 


We  previously  [4]  looked  at  the  properties  of  the  matrix  K  and  its  singular  values  for  a  medium  with  non-zero  attenuation. 
We  found  that  2l(q,  z)  was  much  more  reliably  recovered  for  q  with  magnitude  smaller  than  the  wave  number  -  than  for 
larger  q  .  We  also  observed  that  smaller  eigenvalues  typically  corresponded  to  eigenvectors  with  most  of  their  energy  at 
greater  depths.  The  behavior  of  eigenvectors  and  eigenvalues  indicated  that  shallower  objects  can  be  recovered  much  better 
than  deeper  objects  because  they  correspond  to  eigenvectors  with  much  higher  eigenvalues. 

2.2.  3D  Simulated  pressure  data 

We  simulated  the  3D  attenuated  pressure  data  in  planar  geometry  using  two  methods.  The  SVD-based  method  directly 
computes  the  attenuated  pressure  signals.  The  Eourier-based  method  first  calculates  the  ideal  pressure  and  then  computes 
the  attenuated  pressure  from  it.  We  used  these  two  methods  to  make  sure  that  the  results  were  consistent. 

2.2.1.  SVD-based  simulated  pressure 

We  used  Eqn.  (7)  to  calculate  the  attenuated  pressure.  This  was  implemented  as  follows: 


1.  Take  the  2D  spatial  Eourier  transform  of  the  optical  absorption  function  in  the  x  and  y  directions  to  obtain  A{q,  z). 

2.  Construct  the  matrix  iT„(q,  z)  as  a  function  of  discrete  w,  q  and  z  based  on  Eqn.  (9). 

3.  Multiply  A(q,  z)  by  iT„(q,  z)  and  sum  over  z  to  obtain  P„(q). 

4.  Take  the  inverse  Eourier  transform  of  P„(q)  to  obtain  attenuated  pressure  p(r,  t). 

2.2.2.  Fourier-based  simulated  pressure 

We  used  the  Eourier-based  expression  described  by  Kostli  et  al.  [10]  to  calculate  the  3D  ideal  pressure.  This  expression 
relates  the  3D  Eourier  transform  of  the  optical  absorption  function  to  the  pressure  as: 


1 

(27r)3 


(fkA{k)  cos(a;f)  exp(*k-r), 


(12) 


where,  uj  =  cq  y  -\-  ky-\-  We  evaluated  the  ideal  pressure  on  the  plane  z  =  0.  We  then  used  the  integral  expression 
derived  by  La  Riviere  et  al.  [3]  to  compute  the  attenuated  pressure  from  the  ideal  pressure  data.  Their  expression  is  given 
by: 


p(r,w)  =  /(w) 


cq 

-(w) 


iaocosffn(w) 


Pideaiijc,t)exp^ 


Cq 

Alo) 


■  taoCo|w| 


dt. 


(13) 


2.3.  Image  Reconstruction  Details 

The  images  were  reconstructed  using  the  SVD-based  expression  using  the  two  different  simulated  pressures.  This  was  done 
by  performing  the  SVD  of  matrix  K  to  obtain  the  regularized  pseudoinverse  of  matrix  M.  The  matrix  K  was  zero-padded 
in  the  temporal  frequency  domain.The  mean  value  of  the  eigenvalues  of  M  for  a  specific  q  was  used  to  regularize  the 
pseudoinverse  of  M  for  that  q.  The  images  were  reconstructed  on  a  128x128x32  grid.  The  number  of  transducer  positions 
in  the  X  —  Y  plane  were  128x128  and  32  temporal  samples  were  used.  The  pixel  size  in  the  X  —  Y  plane  was  set  to  0.75 
mm.  The  2-slice  thickness  was  0.75  mm.  The  maximum  temporal  frequency  was  set  to  1.0  MHz.  The  temporal  sampling 
interval  was  set  to  500  nsec.  The  attenuation  coefficient  was  set  to  cm  'rad  's  .  The  speed  of  sound  cq  was  set 

to  1500  m/s  at  =  1  MHz.  A  point  source  was  placed  at  different  depths  and  its  image  was  reconstructed  to  obtain  depth 
resolution. 


3.  RESULTS 

3.1.  Noiseless  Data 

The  reconstructed  images  using  SVD-based  reconstruction  are  shown  below  in  figure  1 .  The  reconstructed  images  bear  a 
close  resemblance  to  the  actual  phantom.  There  are  some  spurious  artifacts  in  the  reconstructed  images  due  to  inter-planar 
interference  and  the  inexact  nature  of  reconstruction  due  to  regularization  of  the  pseudoinverse  of  matrix  M. 


Figure  1.  Reconstructed  2-slices  (0.15  cm,  0.225  cm,  0.3  cm  and  0.45  cm)  of  a  3D  phantom  using  SVD-based  image  reconstruction 
approach  top:  Fourier-based  attenuated  pressure,  middle:  SVD-based  attenuated  pressure,  bottom:  actual  phantom 


3.2.  Noisy  Data 

Gaussian  noise  with  zero  mean  with  10%  amplitude  was  added  to  the  pressure  data  to  obtain  noisy  data.  The  following 
images  were  reconstructed  using  a  bigger  transducer  size  of  200x200  and  a  maximum  frequency  of  10  MHz  since  we 
wanted  to  evaluate  the  performance  of  the  SVD-based  algorithm  at  greater  depths.  Figure  2  shows  the  reconstructed 
phantom  using  noiseless  and  noisy  attenuated  pressure. 


(a)  True  (b)  Noiseless  (c)  Noisy 

Figure  2.  Select  reconstructed  z-plane  at  a  depth  of  1.6875  cm  reconstructed  using  noiseless  and  noisy  SVD-based  attenuated  pressure 


Figure  3  shows  profiles  corresponding  to  the  line  y  =  1.8375  cm  in  the  slice  z  =  1.6875  cm  of  the  reconstructed  A(r). 
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Figure  3.  Line  profile  through  phantom  at  depth  z  =  1.6875  cm  and  y  =  1.8375  cm  for  noiseless  and  noisy  data 

The  reconstructed  images  based  on  the  noisy  and  noiseless  simulated  pressures  obtained  via  SVD-based  method  are  in 
good  agreement  with  the  true  phantom.  The  SVD-based  image  reconstruction  algorithm  offers  comparable  images  even 
with  noisy  data  at  greater  depths. 


3.3.  Depth  Resolution 

Figure  4  shows  the  profile  through  a  point  source  placed  at  various  depths  reconstructed  using  SVD-based  attenuated 
pressure.  The  SVD-based  image  reconstruction  algorithm  shows  very  good  depth  resolution. 
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Figure  4.  Line  profile  of  a  point  source  placed  at  various  depths 


4.  CONCLUSIONS 

We  were  able  to  use  a  previously  derived  operator  relating  the  attenuated  pressure  to  optical  absorption  function  to  recon¬ 
struct  images  in  PAT  in  an  attenuating  medium  for  planar  geometry.  The  image  reconstruction  algorithm  based  on  the  S  VD 
of  this  operator  shows  very  good  depth  resolution  and  noise  stability.  Further  studies  will  include  using  this  algorithm  on 
actual  attenuated  pressure  data. 
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Photoacoustic  image  reconstruction  in  an  attenuating 
medium  using  singular-value  decomposition 

Dimple  Modgil  and  Patrick  J.  La  Riviere 


I.  Introduction 

Photoacoustic  tomography  (PAT)  or  optoacoustic  tomogra¬ 
phy  reconstructs  an  image  of  the  optical  absorption  function 
of  a  medium  using  the  ultrasound  signals  induced  when  it 
is  exposed  to  pulsed  electromagnetic  radiation.  These  signals 
are  generated  due  to  the  thermal  expansion  caused  by  local¬ 
ized  heating  from  the  absorption  of  pulsed  electromagnetic 
radiation  [1].  The  majority  of  the  image  reconstruction  algo¬ 
rithms  in  PAT  so  far  have  assumed  a  non-dispersive  acoustic 
medium.  The  effect  of  frequency-dependent  attenuation  on 
acoustic  waves  can  be  significant  since  PAT  uses  broadband 
detection.  Reconstructed  images  may  exhibit  distortion  and 
artifacts  if  these  effects  are  not  taken  into  account.  Previous 
work  on  dispersive  acoustic  media  done  by  our  group  [2] 
focused  on  incorporating  the  frequency-dependent  attenuation 
effects  into  the  PAT  model.  In  this  paper,  we  will  use  an 
approach  similar  to  that  by  Schotland  et  al.  [3]  [4],  in  optical 
diffusion  tomography,  to  derive  an  inversion  formula  for  the 
optical  absorption  function  using  singular-value  decomposition 
(SVD).  This  formula  is  applicable  in  a  planar  geometry  where 
the  array  of  transducers  lies  in  a  plane.  It  provides  an  insight 
into  the  conditioning  of  the  inverse  problem  and  offers  a 
promising  method  for  image  reconstruction  in  an  attenuating 
medium. 


II.  Methods 

Under  the  constraints  of  negligible  thermal  diffusion,  the 
photoacoustic  pressure,  p(r,  t)  satisfies  the  wave  equation  [5] 
[6]: 


where  cq  is  the  speed  of  sound  in  the  medium,  Cp  is  the 
specific  heat,  and  (3  is  the  coefficient  of  thermal  expansion. 
H{r,t)  is  the  heating  function  that  denotes  the  energy  de¬ 
posited  per  unit  time  per  unit  volume  in  the  medium  by  the 
illuminating  optical  pulse.  One  can  express  H{r,t)  as: 


vMr,w)  +  ^p(r,w)  =  *^A(r)/(w),  (2) 

Cq  Op 

where  uj  is  the  temporal  frequency. 

In  the  presence  of  attenuation,  the  wave  number  becomes 
complex  and  is  given  by  [7]: 

k{uj)  =  +ia{u;).  (3) 

c(w) 

For  tissues,  the  ultrasonic  attenuation  a  is  given  to  a  good 
approximation  by: 


a{Lo)  =  ao\uj\.  (4) 

Using  Kramers-Kronig  relations,  one  obtains  the  depen¬ 
dence  of  phase  velocity  c{ui)  on  frequency  due  to  frequency- 
dependent  attenuation  as: 


1 

c(w) 


1 

Co 


-aoln\  —  |. 
TT  Wq 


So  for  an  attenuating  medium,  eqn.  (2)  becomes: 


(5) 


\/^p{r,uj)  +  k'^{uj)p{r,u))  =  - A{r)  I  (u) .  (6) 

Op 

This  equation  can  be  solved  using  standard  Green’s  function 
methods  as  [6]: 


p{r,uj)  =  —iriujl{uj)  J  J  J  d^i"'A{r')G{r  —  r'),  (7) 

where  r)  = 

A.  Angular  spectrum  expansion  of  the  measured  pressure 
signals 

The  angular  spectrum  expansion  of  the  Green’s  function  is 
given  by  [8]: 


H{r,t)  =  A{r)I{t), 

where  A(r)  is  the  optical  energy  absorption  at  r  and  I{t) 
is  the  temporal  profile  of  the  optical  pulse.  On  taking  the 
Fourier  transform  of  the  wave  equation  with  respect  to  time, 
one  obtains: 
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G(r)  = 


exp{iKr) 

47rr 
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X  exp  l^l  ^of  +  ci^  —  kf  +  ixux  +  iyo-^  daxduy. 

(8) 


Substitute  this  in  eqn.  (7)  and  consider  the  pressure  measure¬ 
ments  on  the  plane  z  =  0  and  we  get: 


p{x,y,0,uj)  = 


—iriujl{uj) 

Stt^ 


X  I  I  da^day 


d^r'A{x',y',z') 
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^al  +  al-k"^ 

X  exp  (^—\z'\^a'l  +  a^  —  kP'  +  ix' a ^  +  iy' . 


(9) 


We  assume  that  the  photoacoustic  object  lies  in  the  plane 
z'  >  0.  On  taking  the  Fourier  transform  on  both  sides  and 
reducing  the  resulting  expression  one  obtains  the  angular 
spectrum  expansion  of  measured  photoacoustic  pressure  as: 


p{kx,ky,uj)=  - ,  /  dz'A{kx,ky,z') 

X  exp  (^-z'^kl  +  fcy  -  k?‘{ijj)^  .  (10) 

1)  SVD  of  integral  operator  relating  pressure  to  optical 
absorption  coefficient:  We  follow  Schotland’s  [3]  [4]  approach 
to  optical  diffusion  tomography  to  obtain  an  integral  operator 
relating  the  measured  attenuated  pressure  to  the  optical  ab¬ 
sorption  coefficient. 

Define  two-dimensional  spatial  wave  vector  q  as: 

q={k^,ky,0).  (11) 

Discretize  the  temporal  frequency,  uj  as  w„.  One  can  then 
write  the  pressure  as: 


p(q,w„) = 


-iriujni  {u)„) 


dzA{q,  z) 


-  fc2(w„) 

X  exp  [-zsjq^  -  kffuJrS)  ■ 
Define  P„(q)  =  p(q, a;„)  ,  kn  =  k{uin)  and 

Thus,  equation  (12)  becomes: 


(12) 


(13) 


-Pn(q)  =  J  dzA{q,z)Kn{q,z).  (14) 

This  equation  can  be  inverted  as  [4]: 

A{q,z)  =  (15) 

m^n 

where 

MmniA}  =  J  dzKjn{q,z)K*{q,z).  (16) 

Taking  the  inverse  Fourier  transform  of  eqn.  (15)  with 
respect  to  q  gives: 

\  m,n 

(17) 


Define  (5„(q)  =  This  function  (5„(q)  is  com¬ 

plex  for  non-zero  attenuation.  One  can  simplify  equation  (16) 
to  obtain: 


Jo  ^  (87r2)2Q„(q)Q*(q) 

X  exp  (-z(Qm(q)  +  (5^(q))) 


This  can  be  reduced  to: 


(18) 


AT  /  \  V  ujjjiUjjiI fuj  m  )n^n)  (  1  \ 

(87r2)2g„(q)g*  (q)  (Q^(q)  +  Q*  (q))  j  ' 

(19) 

The  computation  of  the  inverse  of  this  matrix  is  the  key  step 
in  the  procedure  to  recover  the  optical  absorption  function.  Its 
pseudoinverse  can  be  computed  by  performing  SVD  of  matrix 
K.  Some  of  the  eigenvalues  of  this  matrix  may  be  zero  or  very 
small.  One  will  need  to  use  regularization  to  circumvent  this 
problem. 

Thus,  the  procedure  for  recovering  xl(r)  is: 

1)  Take  the  Fourier  transform  of  the  measured  pressure 
(at  z  =  0)  with  respect  to  time  and  the  2D  spatial 
components  to  obtain  p{kx,ky,u}). 

2)  Discretize  asp(q,  w„). 

3)  Compute  the  matrix  Kn{q,z)  using  eqn.  (13)  . 

4)  Compute  the  pseudoinverse  of  matrix  Mmn{q)  by  per¬ 
forming  SVD  of  matrix  K. 

5)  Compute  ^(q,  z)  using  eqn.  (15). 

6)  Take  the  2D  inverse  Fourier  transform  of  Tl(q,  z)  to 
obtain  A(r). 

2 )  Calculation  of  eigenvalues  of  M  matrix  with  non-zero 
attenuation:  We  looked  at  the  properties  of  the  matrix  K  and 
its  singular  values  for  a  medium  with  non-zero  attenuation. 
We  computed  this  matrix  for  a  tissue  with  attenuation  ag  = 
i2;^cm‘'rad‘'s.  A  set  of  32  discrete  temporal  frequencies  and 
64  discrete  spatial  frequencies  was  used.  Temporal  frequencies 
ranged  from  0  to  12.17  MHz.  and  spatial  frequencies  ranged 
from  0  to  82.23  cm  '.  We  performed  an  SVD  of  the  matrix 
K.  The  eigenvalues  of  the  matrix  M  were  obtained  from  the 
singular  values  of  matrix  K. 

The  surface  plot  of  eigenvalues  of  M  vs.  the  spatial  fre¬ 
quency  q  is  shown  in  figure  1.  From  this  plot,  one  can  see  that 
values  of  2D  spatial  frequency  vector  q  that  are  much  smaller 
than  the  wave  vector  are  recovered  much  better  than  the  ones 
that  are  of  higher  magnitude.  Figure  2  shows  the  variation 
of  eigenvectors  of  K^K  with  z.  In  this  figure,  the  vertical 
direction  is  increasing  depth  and  the  horizontal  direction 
is  decreasing  eigenvalues.  Figure  3  shows  the  variation  of 
the  eigenvectors  of  matrix  K^K  with  depth  for  a  specific 
eigenvalue.  In  these  plots,  a  lower  column  number  corresponds 
to  a  higher  eigenvalue.  From  these  plots  one  concludes  that: 

1)  The  values  of  2D  spatial  frequency  wave  vector  ’q’  that 
are  smaller  than  the  wave  vector  are  recovered  much 
better. 

2)  The  eigenvectors  that  are  non-zero  at  greater  depths 
correspond  to  smaller  eigenvalues. 


The  behavior  of  eigenvectors  and  eigenvalues  indicates  that 
shallower  objects  can  be  recovered  much  better  than  deeper 
objects  because  they  correspond  to  eigenvectors  with  much 
higher  eigenvalues. 
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Figure  1.  Plot  of  eigenvalues  of  M  vs.  q 


Figure  2.  Real-part  of  a  set  of  eigenvectors  of  matrix  K  for  a  specific 
q,  x-axis:  decreasing  eigenvalues,  j/-axis:increasing  depth 


B.  3D  Simulated  pressure  data 

We  simulated  the  3D  attenuated  pressure  data  in  planar 
geometry  using  two  methods.  The  SVD-based  method  directly 
computes  the  attenuated  pressure  signals.  The  Fourier-based 
method  first  calculates  the  ideal  pressure  and  then  computes 


(a)  Column=10 


(c)  Column=230 

Figure  3.  Variation  of  select  eigenvectors  with  depth 


the  attenuated  pressure  from  it.  We  used  these  two  methods 
to  make  sure  that  the  results  were  consistent  and  to  avoid  the 
possibility  of  an  inverse  crime. 

1)  SVD-based  simulated  pressure:  We  used  eqn.  (14)  to 
calculate  the  attenuated  pressure.  This  was  implemented  as 
follows: 

1)  Take  the  2D  spatial  Fourier  transform  of  the  optical 
absorption  function  in  the  x  and  y  directions  to  obtain 
i(q,2). 

2)  Construct  the  matrix  iT„(q,  z)  as  a  function  of  discrete 
w,  q  and  2  based  on  eqn.  (13). 

3)  Multiply  2l(q,  z)  by  Ff„(q,  z)  and  sum  over  z  to  obtain 
-Pn(q)- 

4)  Take  the  inverse  Fourier  transform  of  P„(q)  to  obtain 
attenuated  pressure  p{r,t). 

2 )  Fourier-based  simulated  pressure:  We  used  the  Fourier- 
based  expression  described  by  Kostli  et  al.  [9]  to  calculate 
the  3D  ideal  pressure.  This  expression  relates  the  3D  Fourier 
transform  of  the  optical  absorption  function  to  the  pressure  as: 


Pi'T,*) 


cos(a;t)  exp(ik-r),  (20) 


(27r)3 

where,  w  =  +  fcf .  We  evaluated  the  ideal  pres¬ 

sure  on  the  plane  z  =  0.  We  then  used  the  integral  expression 
derived  by  La  Riviere  et  al.  [2]  to  compute  the  attenuated 
pressure  from  the  ideal  pressure  data.  Their  expression  is  given 
by: 


p{r,uj)  =  i{uj) 


Co 
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(21) 


C.  Image  Reconstruction  Details 

The  images  were  reconstructed  using  the  SVD-based  ex¬ 
pression  using  the  two  different  simulated  pressures.  This 
was  done  by  performing  the  SVD  of  matrix  K  to  obtain  the 
regularized  pseudoinverse  of  matrix  M.  The  mean  value  of  the 
eigenvalues  of  M  for  a  specific  q  was  used  to  regularize  the 
pseudoinverse  of  M  for  that  q.  The  images  were  reconstructed 
on  a  128x128x32  grid.  The  number  of  transducer  positions  in 
the  X  —  Y  plane  were  128x128  and  32  temporal  samples  were 
used.  The  pixel  size  in  the  X  —  Y  plane  was  set  to  0.75  mm. 
The  z-slice  thickness  was  0.75  mm.  The  maximum  temporal 
frequency  was  set  to  1.0  MHz.  The  temporal  sampling  interval 
was  set  to  500  nsec.  The  attenuation  coefficient  was  set  to 
ao  =  cm‘'rad‘'s  .  The  speed  of  sound  cg  was  set  to 

1500  m/s’^at  =  1  MHz. 


HI.  Results 

The  reconstructed  images  using  SVD-based  reconstruction 
are  shown  below  in  figure  4.  The  reconstructed  images  bear 
a  close  resemblance  to  the  actual  phantom.  There  are  some 
spurious  artifacts  in  the  reconstructed  images  due  to  inter- 
planar  interference  and  the  inexact  nature  of  reconstruction 
due  to  regularization  of  the  pseudoinverse  of  matrix  M. 


Figure  4.  Reconstructed  2:-slices  (0.15  cm,  0.225  cm  and  0.30  cm)  of  a  3D 
phantom  using  SVD-based  image  reconstruction  approach  top:  Fourier-based 
attenuated  pressure,  middle:  SVD-based  attenuated  pressure,  bottom:  actual 
phantom 


True 

SVD-based  pressure 
Fourier-based  pressure 


Figure  5.  Profile  along  the  line  y  =  4.725  cm,  2:  =  0.225  cm,  solid  line 
-  true  profile,  dotted  line  -  using  SVD-based  pressure,  dashed  line  -  using 
Fourier-based  pressure 


Figure  5  shows  profiles  corresponding  the  line  y  =  4.725 
cm  in  the  slice  z  =  0.225  cm  of  the  reconstructed  7l(r). 

The  reconstructed  images  based  on  the  simulated  pressures 
obtained  via  SVD-based  and  Fourier-based  methods  are  in 
good  agreement  with  the  true  phantom. 

IV.  Discussion 

We  have  derived  an  SVD-based  algorithm  to  reconstruct  the 
optical  absorption  function  in  PAT  in  an  attenuating  medium. 
We  looked  at  the  eigenvalues  of  the  matrix  M  which  is  used 
to  recover  the  optical  absorption  function.  We  found  that 
smaller  values  of  2D  spatial  frequency  vector  q  are  recovered 
much  better  than  the  larger  values.  The  SVD  approach  to 


the  attenuation  problem  in  PAT  offers  a  promising  method 
for  image  reconstruction  that  is  direct  and  computationally 
efficient.  It  also  provides  a  way  to  study  the  conditioning  of 
this  inverse  problem.  Further  studies  need  to  be  conducted 
with  actual  attenuated  pressure  data.  We  also  need  to  perform 
noise  and  resolution  studies  to  evaluate  the  image  quality  of 
images  reconstructed  using  the  SVD-based  image  reconstruc¬ 
tion  method. 
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